act_site = [1 2 4 7 8];
pas_site = [3 5 6 9];
con_site = 10;

for site = 1:9
    for layer = 1:4
        for time = 1:4
            cur_fb_act = proj_meta(site).rd(layer,time).act(:,1:7500);
            cur_pb_act = proj_meta(site).rd(layer,time).act(:,7501:15000);

            % generate binary activity during feedback.
            cur_fb_bi = proj_meta(site).rd(layer,time).velM_smoothed(1:7500);
            cur_fb_bi = -cur_fb_bi; % make the velocity of forward-running positive.
            cur_fb_bi(cur_fb_bi > 0.001) = 1;
            cur_fb_bi(cur_fb_bi <= 0.001) = 0;
            
            % generate binary activity during playback.
            cur_pb_bi = proj_meta(site).rd(layer,time).velM_smoothed(7501:15000);
            cur_pb_bi = -cur_pb_bi; % make the velocity of forward-running positive.
            cur_pb_bi(cur_pb_bi > 0.001) = 1;
            cur_pb_bi(cur_pb_bi <= 0.001) = 0;
            
            % generate the frame id of ps onset.
            cur_ps_ons = find(diff(proj_meta(site).rd(layer,time).ps_id(1:7500))>1)+1;
            % remove early & late ps.
            cur_ps_ons(cur_ps_ons > 7425) = [];
            cur_ps_ons(cur_ps_ons < 6) = [];
            
            % generate the frame id of playback halt onset.
            cur_pb_ons = find(diff(proj_meta(site).rd(layer,time).ps_id(7501:15000))>1)+1;
            % remove early & late pb.
            cur_pb_ons(cur_pb_ons > 7425) = [];
            cur_pb_ons(cur_pb_ons < 8) = [];
            
            rps_nbr = 0; % number of real ps
            rps_id = [];
            for ps_id = 1:length(cur_ps_ons)
                if sum(cur_fb_bi(cur_ps_ons(ps_id)-5:cur_ps_ons(ps_id)+15)) == 21
                    rps_nbr = rps_nbr+1;
                    rps_id(rps_nbr) = ps_id;
                    all_rps_nbr(site,time) = rps_nbr;
                end
            end
            
            rpb_nbr = 0; % number of real pb halt
            rpb_id = [];
            for pb_id = rps_id % rpb is selected based on rps
                if sum(cur_pb_bi(cur_pb_ons(pb_id)-7:cur_pb_ons(pb_id)+22)) == 0
                    rpb_nbr = rpb_nbr+1;
                    rpb_id(rpb_nbr) = pb_id;
                    all_rpb_nbr(site,time) = rpb_nbr;
                end
            end
            
%             if rps_nbr > 10
%                 rps_id = rps_id(1:10);
%             end
%             
%             if rpb_nbr > 10
%                 rpb_id = rpb_id(1:10);
%             end
            
            for cell_id = 1:size(cur_fb_act,1)
                rps_resp_cnt = 0;
                rps_resp = [];
                for rps_act_nbr = 1:length(rps_id)
                    rps_frame_id = cur_ps_ons(rps_id(rps_act_nbr));
                    rps_resp_cnt = rps_resp_cnt+1;
                    rps_resp(rps_resp_cnt,:) = cur_fb_act(cell_id,rps_frame_id-7:rps_frame_id+75);
                    cur_peak_fb_act = cur_fb_act(cell_id,rps_frame_id+16:rps_frame_id+30);
%                     cur_pre_fb_act = 0;
                    cur_pre_fb_act = cur_fb_act(cell_id,rps_frame_id-7:rps_frame_id);
                    prps_resp{site,layer,time,cell_id,rps_resp_cnt} = mean(cur_peak_fb_act)-mean(cur_pre_fb_act);
                end
                rps_resp_by_cell{site,layer,time,cell_id} = rps_resp;
            end
            
            for cell_id = 1:size(cur_pb_act,1)
                rpb_resp_cnt = 0;
                rpb_resp = [];
                for rpb_act_nbr = 1:length(rpb_id)
                    rpb_frame_id = cur_pb_ons(rpb_id(rpb_act_nbr));
                    rpb_resp_cnt = rpb_resp_cnt+1;
                    rpb_resp(rpb_resp_cnt,:) = cur_pb_act(cell_id,rpb_frame_id-7:rpb_frame_id+75);
                    cur_peak_pb_act = cur_pb_act(cell_id,rpb_frame_id+16:rpb_frame_id+30);
%                     cur_pre_pb_act = 0;
                    cur_pre_pb_act = cur_pb_act(cell_id,rpb_frame_id-7:rpb_frame_id);
                    prpb_resp{site,layer,time,cell_id,rpb_resp_cnt} = mean(cur_peak_pb_act)-mean(cur_pre_pb_act);
                end
                rpb_resp_by_cell{site,layer,time,cell_id} = rpb_resp;
            end
        end
        cell_cnt{site,layer} = size(proj_meta(site).rd(layer,1).act,1);
    end
end

for time = 1:4
    for rps_nbr = 1:33
        rps_act{time,rps_nbr} = mean([prps_resp{act_site,:,time,:,rps_nbr}]);
        rps_pas{time,rps_nbr} = mean([prps_resp{pas_site,:,time,:,rps_nbr}]);
    end
    
    for rpb_nbr = 1:22
        rpb_act{time,rpb_nbr} = mean([prpb_resp{act_site,:,time,:,rpb_nbr}]);
        rpb_pas{time,rpb_nbr} = mean([prpb_resp{pas_site,:,time,:,rpb_nbr}]);
    end
end

for time = 1:4
    rps_act_seq{time} = horzcat(rps_act{time,:});
    rps_pas_seq{time} = horzcat(rps_pas{time,:});
end

for time = 1:4
    rpb_act_seq{time} = horzcat(rpb_act{time,:});
    rpb_pas_seq{time} = horzcat(rpb_pas{time,:});
end

figure;hold on
plot(smooth2(horzcat(rps_act_seq{:}),3))
plot(smooth2(horzcat(rps_pas_seq{:}),3),'r')

figure;hold on
plot(smooth2(horzcat(rpb_act_seq{:}),3),'--')
plot(smooth2(horzcat(rpb_pas_seq{:}),3),'r--')

for rps_nbr = 1:33
    rps_act_all(rps_nbr) = mean(vertcat(rps_act{:,rps_nbr}));
    rps_pas_all(rps_nbr) = mean(vertcat(rps_pas{:,rps_nbr}));
end

for rpb_nbr = 1:22
    rpb_act_all(rpb_nbr) = mean(vertcat(rpb_act{:,rpb_nbr}));
    rpb_pas_all(rpb_nbr) = mean(vertcat(rpb_pas{:,rpb_nbr}));
end

figure;hold on
plot(smooth2(rps_act_all,3))
plot(smooth2(rps_pas_all,3),'r')
plot(smooth2(rpb_act_all,3),'--')
plot(smooth2(rpb_pas_all,3),'r--')